The Asymmetric Superfluid Local Density Approximation (ASLDA) 
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Here we describe the form of the Asymmetric Superfluid Local Density Approximation (aslda), 
a Density Functional Theory (dft) used to model the two-component unitary Fermi gas. We give 
the rational behind the functional, and describe explicitly how we determine the form of the dft 
from the to the available numerical and experimental data. 



I. INTRODUCTION 

Here we describe the formulation of the Asymmet- 
ric Superfluid Local Density Approximation (aslda), 
which is a Density Functional Theory (dft) describing 
normal and superfluid systems comprising two species 
of fermion. We show how the Monte-Carlo data of 
Refs. 0, [2, [3, i, S, H, 0, H was incorporated into the 
ASLDA dft used in the papers d, [lij ■ In its present for- 
mulation, the ASLDA functional describes two species of 
fermions - denoted a (spin-up) and b (spin-down) - with 
equal masses m a — rtib — rn, and interacting through a 
resonant attractive inter-species s-wave two-body inter- 
action described solely by the infinite two-body scatter- 
ing length a s = oo. Throughout this paper, we use units 
where H = m = 1 to simplify the notation. 

The model is tuned to reproduce the thermodynamics 
of the homogeneous normal phase, and the homogeneous 
symmetric (n a — rn,) superfluid phase (sf). Properties 
of these phases have been evaluated using accurate non- 
perturbative Monte-Carlo calculations, including the en- 
ergy of the interacting normal state with varying degrees 
of polarization n a ^ rib (Fig. CO, the energy of the SF 
phase £ = £/£ F g = 0.40(1) [1, [|, the quasi-particle 
dispersion relationship in the SF phase (Fig. [2]) [2| , and 
the effective mass of a spin-down fermion immersed in 
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x = Ub/n a 

FIG. 1: (Color online) The dimensionless convex func- 
tion g(x) [IIj that defines the energy density £(n a ,nb) = 

l^r( 67I " 2 ) 2/3 K5 , ( a; )] 5/3 - The points with error-bars (blue 
online) are the Monte-Carlo data from Refs. @, [1 @|. The 
fully-paired solution g(l) = (2£) 3//5 is indicated to the bot- 
tom right, and the recent MIT data [l^] is shown (light x) for 
comparison. The phase separation discussed in Refs. [U, ID, |(| 
is shown by the Maxwell construction (thin black dashed line). 



a sea of spin- up fermions 0, [H, Q. We describe here 
how to incorporate this high-quality information into the 
ASLDA functional, allowing it to accurately describe all 
the homogeneous normal and superfluid properties of the 
unitary Fermi gas. The ASLDA functional thus provides 
means of using the quantitative non-perturbative infor- 
mation about homogeneous phases to explore the struc- 
ture of inhomogeneous systems. This has led to the pre- 
diction of a non-trivial supersolid Larkin-Ovchinikov (lo) 
phase [10( of the polarized unitary Fermi gas. 



II. THE FUNCTIONAL 

We choose to use an ASLDA functional with the follow- 
ing form, 



d 3 r <^ £[n atb {r),T a ,b(r),v(r)] + 



+ V cxt (r)[n a (r) + n b (r)] + 



sources . 



(1) 



where the (local) energy-density £ is a function the 
densities. We construct the densities from the two- 
component Bogoliubov (Nambu-Gor'kov) quasiparticle 




FIG. 2: Fit of the Monte-Carlo data for the quasiparticle dis- 
persions from [2| with the BCS form ([5]). This, along with £ 
is used to determine the parameters of the fully paired super- 
fluid phase SF. 
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wavefunctions that satisfy the generalized Kohn-Sham 
equation Hks ' VVi = E n ip n (see below): 



u n (r) 
v n (r) 



(2) 



The densities are constructed from these by summing 
over all eigenstates E n weighted with the Fermi distri- 
bution function fp(E), where j3 = 1/T. (A small finite 
temperature can help the initial part of the self-consistent 
calculations to converge, see [13].) 



n a (r) 


=E 

n 


u„(r)| 2 //3(£«), 


(3a) 


n b (r) 


=E 

n 


v n (r)\ 2 M-E n ), 


(3b) 


T a(r) 


=E 

n 


Vu n (r)\ 2 f p (E n ), 


(3c) 


n(r) 


=E 


Vv n {r)\ 2 f p {-E n ), 


(3d) 


v{r) 


=E 




(3e) 



We express our functional in terms of these densities 
(from here on we drop the explicit r dependence to sim- 
plify the notation) 

£[n a ,b, i~ a ,b, v\ = 

= a a (n a ,n b )Y + a b (n a ,n b )^+ 

D(n a , rib) + g e s{n a ,n b )v*v (4a) 

where a a ^ b {\n a ,Xn b ) — a(n a ,n b ) are the inverse ef- 
fective masses, D(\n a ,\n b ) — \ 5 ^ 3 D(n a ,n b ) contains 
the density dependent portion of the functional and 
C(An a , \n b ) — A 1 / 3 C(n a , n b ) is an appropriate density 
dependent inverse effective coupling, with regulator A: 



1 

ffcff 



= C(n a ,n b ) - A. 



(4b) 



The divergences appear in both the pairing and kinetic 
terms such that A = —g c s^ and the combination 



a b - 



n 



(4c) 



are finite. Thus, the energy functional may be expressed 
in terms of finite combinations 



a b - 



n 



D 



(4d) 



where A*i> = i>*A — —g e si/*i/. One recovers a typical 
Kohn-Sham functional by setting a ajb = 1 and A = 
(y = 0). Regularization is required to describe super- 
fluids (A ^ and v ^ 0) because the anomalous den- 
sity v and the kinetic energy densities r a b diverge. This 
regulator dependence is removed when one considers the 



proper combinations of terms described above. Allow- 
ing the effective mass to deviate a ^ 1 in the SLDA de- 
scription of the trapped fermions improved the agreement 
with Monte-Carlo (mc) results [|, til, EE 03 . and can 
be constrained by non-perturbative results, so we keep 
this generalization. The coefficient functions a a {n a ,n b ), 
a b (n a , rib), C(n a , n b ), and D{n a , n b ) are chosen to match 
the thermodynamic properties of the system in the ho- 
mogeneous limit. The lack of scales in the unitary limit 
further constrains the forms of these coefficients to be 
homogeneous functions of specific degree in the densities 
n a ,b- 

We established this form of the ASLDA functional using 
the following guiding principles: 

Gradient Expansion: We start by formulating a func- 
tional applicable to slowly varying systems. To this 
end, we have neglected gradient terms. Indeed, the 
symmetric SLDA functional has had remarkable suc- 
cess without any gradient corrections Q , so we as- 
sume that this property holds for asymmetric sys- 
tems and start with only local terms. An important 
future direction will be to quantify the effects of 
gradient corrections, and to extract the coefficients 
of the leading gradient terms. 

Note that the success of the SLDA implies that the 
subset of gradient corrections included implicitly 
through the standard kinetic terms provides the 
dominant gradient contribution. These corrections 
are also included in the ASLDA. 

To justify the omission of gradient terms (beyond 
those contained in T a ,b)i consider the lowest order 
gradient correction to the SLDA functional: this is 
of the form |Vn| 2 /n and would lead to a correc- 
tion in the total energy of a system of N har- 
monically trapped fermions that scales as N 2 / 3 . 
The extremely good agreement between the SLDA 
functional in Ref. and the ab initio results of 
Refs. [H, [IE EH indicate that the coefficient of 
this correction must be extremely small. One can 
also show that, in the dilute limit, the strength of 
the gradient terms is controlled by the size of the 
length scales such as the effective range and p-wave 
scattering length of the interaction [13, EB] • In the 
unitary Fermi gas, both of these vanish. 

In addition, the leading gradient corrections that 
could appear away from the symmetric limit x = 
rib/n a = 1 are suppressed by (I — x) 2 . Thus, we 
expect them to have little effect on phases close to 
the fully paired SF state, for example, the LO phase 
discussed in (foj . 

Finally, Galilean invariance requires that if a a , b 1 
then r a ,b be replaced by Tah-p 2 b /n a ^ b where p a _ b is 
the local current density . This correction is only 
required to discuss states that break time-reversal 
invariance (i.e. that contain currents). 

Simplicity: Due to the fermion sign problem, there are 
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very few reliable calculations of properties in the 
polarized regime, thus there is not enough data 
to properly constrain a fully general functional. 
Monte-Carlo simulations, however, have provided 
reasonably constraints on the form of the nor- 
mal state energy density, thus we allow for full 
generality in terms of the functional dependence 
on the densities, but restrict the general depen- 
dence on the anomalous and kinetic densities. The 
most general form of local function would allow 
the unknowns to depend on all of the dimension- 
less combinations of densities, including for exam- 
ple, an arbitrary dependence on the dimensionlcss 
regulator-invariant combination [a a T a /2 + a b r b /2 — 
A*y]/n 5 / 3 . Presently, we see no need for this added 
complication. 

Quasiparticle Dispersion: Monte-Carlo calculations 
about the symmetric phase have suggested that the 
low-temperature quasiparticle dispersions are well 
described by the BCS form 



k 2 



(5) 



where a, (U — fi), and A are effective parameters. 
This form of dispersion follows from introducing the 
anomalous pairing density v through the quadratic 
form v*v along with canonical kinetic terms. 

Decoupling of Paired and Normal States: In our 

actual formulation, we also neglect the general 
density dependence of the pairing interaction 
C(n a ,n b ), replacing this by the same term with a 
single constant used in the SLDA. The justification 
of this is two- fold: 1) The SLDA worked very well, 
even including small polarizations. 2) The energy 
of interaction of the normal state is well approxi- 
mated by considering only the parameters of the 
SLDA without any additional density dependence. 

Together these suggest that the physics of the su- 
pcrfluid state somewhat decouples from the physics 
of the normal state, allowing one to characterize the 
fully paired superfluid independently of the normal 
state. 

This implies a qualitative ansatz of this functional: 
that the structure of the polarized phases arises 
from the competition between the fully paired su- 
perfluid physics and the interacting normal state 
physics. In principle, it is possible that some quali- 
tatively new description is required to properly ac- 
count for the structure of polarized fermionic mat- 
ter - for example to describe the a ppe arance of p- 
wave pairing at large polarizations [19( - but to fix 
such a description will require high precision ab ini- 
tio calculations and/or experiments that have not 
yet been realized. 



Varying this functional with respect to the quasi-particle 
wavefunctions and occupation numbers gives the gener- 
alized Kohn-Sham equation Hks • i>n — E n ip n , where 



H 



KS 



K a - Ma 

A 



U a 



A* 



-Kb + - U b 
and the kinetic and potential operators are 



K a u = --V • {a a Vu) 



U a = 



dc 



dn a 2 
d\na + 



dD 

dn a 



dC 
dn a 



A| 2 + 



+ dn a 
with the notations 

ot a ± a b 
a ± = « 1 



-f ^ V 



C|A| 



(6a) 



(6b) 



(6c) 



T±=T a ±T b , n ± =n a ±n b . (6d) 



The form of these operators for species b are obtained 
by interchanging a «-> b. Note that the terms have been 
grouped so that the ultraviolet (uv) divergences arising 
from the local form of the anomalous density v cancel 
in the last term (see Refs. [t| H(| for details). All other 
terms arc finite. 



III. FITTING THE FUNCTIONAL 

To further specify the functional, we must fix the forms 
of the functions a a fi{n a ,n b ), C(n a ,n b ), and D{n a ,n b ). 
To do this, we characterize the thermodynamic proper- 
ties of the system, which are fortunately quite tightly 
constrained [Tl[, and have both calculational [3. IH ITU l21j 
and experimental fl2l . |22| verification. 

The form of the functional allows us to consider two 
different species (for example, with different masses), but 
we are interested in systems with two identical species. 
The functional must thus exhibit the discrete symmetry 
a «-» b. This constrains the form of the functions a Qj0 , D 
and C: 



a a {n a ,n b ) = a b (n b ,n a ) 
D(n a ,n b ) = D(n bl n a ) 
C(n a ,n b ) = C(n b ,n a ). 



(7a) 
(7b) 
(7c) 



Dimensional analysis determines the overall scaling of the 
functions, and we may fully parametrize the functional 
with three dimensionless functions a(x), b(x), and j(x) 
of the asymmetry x — n b /n a 6 [0, 1], with the comple- 
mentary region determined by symmetries |J7J| : 



a a (n a ,n b ) = a(x), 
a b (n a ,n b ) = a(l/x), 



(8a) 
(8b) 



D{n a ,n b ) = T7r - 5 (3{x), (8c) 



C(n a ,n b ) = 



IOtt 2 
{n a + n b ) 1/3 
7(x) 



(8d) 



4 



As a technical note, we ensure that our parametrization 
is smooth at x = 1 by letting a(x) be smooth over [0, oo), 
and by forming smooth even functions over the variable 
z G [—1, 1] where \ax — ctanh -1 z: thus (3(z) = (3(x(z)) 
and j(z) = j(x(z)) are smooth at x = 1 if and only if 

(3(z) and j(z) are smooth even functions. The choice of 
the parameter c is made so that the interpolations of the 
Monte-Carlo data are well behaved. 

To determine these dimensionless functions, we match 
the functional to the Monte-Carlo calculations of pure 
and homogeneous thermodynamic phases. These phases 
possess no gradients, and so solving the density func- 
tional for homogeneous matter is equivalent to perform- 
ing a simple Thomas Fermi type calculation with the 
added complication that the parameters: m a ^, etc. de- 
pend on the densities, which must be determined self- 
consistently. This gives rise to a set of non-linear set of 
equations that can be fairly easily solved. 

We start with the homogeneous fully paired superfluid 
phase. This is described by the three numbers a = a(l), 
(3 = (3(1), and 7 = 7(1) which may be extracted by fitting 
the quasi-particle dispersion relationship and the energy. 
As discussed earlier, we simplify the functional depen- 
dence of the function 7(2;) = 7 by simply keeping this 
constant. The remaining functional forms a(x) and (3(x) 
are determined by fitting the energy of the homogeneous 
normal state to Monte-Carlo results. This completely 
specifies the functional in a unique manner as we shall 
now describe. 



A. Symmetric Superfluid Properties: 

As suggested in Q , by considering the calculated prop- 
erties of the fully paired symmetric superfluid, one may 
determine the values of the functions a(l), (3(1), and 7(1) 
at the point x = nb/n a = 1. We start by taking 



rib = n — 



(9) 



where n + = n a + rib = 2n is the total density. Di- 
mensional analysis determines the following forms of the 
derivatives (evaluated at n a = rib = n): 



da.4 



dctj 



= 0, 



dn a drib 
dC _dC _ C 
dn a drib in ' 
3D _ 3D _ 5D 
dn a drib 3n 

Thus, the effective potentials U a = Ub 
following simplified form, 



(10a) 
(10b) 
(10c) 

U have the 



U(n, n) 



5D(n,n) C(n,n) 



3n 



■in 



A*A, 



(11) 



and one may take a linear combination of the kinetic 
terms to obtain 



K 4 



K n + K h 



o + V 2 



(12) 



Thus, the symmetric phase depends only on three pa- 
rameters a = a(l), [3 = /3(1), and 7 = 7(1) via: 

a + (n,n) — a(l), (13a) 

1/3 

C(„,„)= + , (13b) 

D(n, n) = (3(1) ^ ^ = (3(1) £ FG . (13c) 

For a given inverse effective mass a, the other two pa- 
rameters 7 and (3 may be fit by requiring that the energy 
and spectral gap satisfy 



£sf = £ (n, n) = ££ FG = £ 

A = r/€F =T) 



(37T 2 71 + ) 5 / 3 

IOtt 2 : 
(37r 2 n+) 2 / 3 



where 



kl (37r 2 n+) 2 / 3 



k F = (37T 2 ?1 + ) 1/3 

£fg 



2 



= 2 



hp 



2 

(3?r 2 ^ 



\5/3 



20tt 2 



IOtt 2 



= -n+e F . 
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(14a) 
(14b) 

(15a) 
(15b) 

(15c) 



The T = symmetric state is characterized by the inte- 
grals 



£sf = 

C = 
where 



d 3 fc 
d 3 fc 



1 - — 

Ek 



1 - 



E k 



2E k 



,V3 



(16a) 

(3£fg, 
(16b) 

(16c) 



(16d) 
(16e) 



Given fixed values of a, £, and 77, we proceed as follows: 
1) Choosing the density n so that ep = 1: This fixes the 
scale and determines A, 2) Use equation (|16a|) to solve 
for the combination U — \x that appears on the right hand 
side through €k/Ek, 3) Use equation (|16b|) and £ to solve 
for (3, and 4) Use equation (|16c|) to solve for 7. 

The parameters £ and 77 = A/ep have been measured 
by several Monte-Carlo techniques [H, H H 0, Hi ■ We take 



r d 3 fc 


1 1 


J (2tt) 3 


2E~ k " o^fc 2 


fc 2 


(U-n), 


ek = a Y + 


Ek=Je 2 k + 


|A| 2 - 
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the following values in our estimates @, [|[ : 
£(n, n) 



£FG(n,n) 



0.40(1), r) 



A 



0.504(24). (17) 



Note that it was incorrectly stated in [9j that a x= \ could 
also be determined through the condition that the chem- 
ical potential satisfy /i = £cf, but a careful examination 
shows that this is implied by (fT4"l) . It is also clear in 
DFT's developed perturbatively [l7|,ll8( that the effective 
mass is arbitrary. In order to fix the effective mass, we 
match the quasi-particle dispersion relationship as deter- 
mined from the Monte-Carlo results 0. The dispersion 
relationships within our density functional have the form 



Ek 



a k 2 
2 kp 



U 



(18) 



Note that the combination U/ep — £ is fixed from (|16a|) 
and depends only on A/ep, so the quasiparticle disper- 
sion relation is sensitive only to A/e^ and the effective 
mass. We add the value (fTT|) as an additional data-point 
in the fit and perform a non- linear least-squares fit. 

The fit to the Carlson-Reddy is shown in Fig. [2] and 
gives the following parameter values: 1 



a x =i 



where 



/m- 1 = 


1.094(17), 


(19a) 


&=l = - 


-0.526(18), 


(19b) 


7" 1 = - 


-0.0907(77) 


(19c) 


U/e F = - 


-0.491(18), 


(19d) 


A/ £f — 


0.493(12), 


(19e) 


a + P = 


0.567(24). 


(19f) 


(3tt 2 ) 2 / 3 
67 




(20) 





v - 



In principle, one should use some form of ab initio cal- 
culation or experimental measurement for polarized sys- 
tems to determine the dependence of the parameters a, 
/?, and 7 on the polarization x = nb/n a . Unfortunately, 
the fermion sign problem has made this difficult and there 
is presently insufficient quality data to perform such a fit. 
Instead, we make the approximation that 



7(2;) = 7(1) = const. 



(21) 



B. Normal State 

The remaining functional forms for a(x) and f3(x) can 
be extracted from properties of the homogeneous nor- 
mal state at T — by explicitly setting the anoma- 
lous density v = 0. This metastable state may be ex- 
plored by choosing a suitable nodal approximation in the 
fixed-node Monte-Carlo (fn-MC) calculations. Such a re- 
striction seems to remove most of the supcrfluid corre- 
lations from the results, but we suspect that one can- 
not completely remove all contamination from the su- 
perfluid state for small polarizations. For this reason, 
we have only included the data for large polarizations 
x = ni,/n a < 0.5 in our fitting of the parameters for the 
functional (see Fig. [1] and Tab. \U\i . 

The energy-density for the normal phase of homoge- 
neous matter has the form 

c . a a (6^ 2 n a ) 5 / 3 a 6 (6^ 2 n fc ) 5 / 3 
o[n a ,n b \ = 1 h D 



20tt 2 20tt 2 
(67r 2 ) 5 / 3 (n a + n b ) 5 / 3 



20tt 2 



-G(x) 



(22) 



where 
G(x) 



(1 + .t) 5 / 3 (1+x 
1 



(1 + x) 5 / 3 



g^(x). 



(23) 



From this relationship, one can uniquely determine the 
functional form for f3(x) given a form for the inverse ef- 
fective mass oe(x), which we shall construct below, and 
the normal state energy density g(x), which has been 
well-constrained by Monte-Carlo data [3| (see Fig. [IJ . 



C. Effective Mass Parametrization: a(x) 

Several different Monte-Carlo calculations have con- 
strained the effective "polaron" mass of a single spin- 
down fermion in a sea of spin-up fermions, m ~ 




1 Here we have performed a simple three-parameter non-linear 
least-squares fit: this has a quality factor Q = 0.52 which is 
quite good. A different analysis would hold A and U fixed to the 
properly determined values, but the method here is consistent 
since the errors are of the same magnitude. 



FIG. 3: Inverse effective mass a = m~g as a function of 
z = tanhMni with c = 1. The function is a parabolic fit 
through the three data-points. 



6 



1.04(3)m US,!!], and the effective mass of the particles 
in the background gas will not be altered in the limit of 
extreme polarization. Thus, the endpoints of the function 
a(x) are constrained: a(0) = 1, a(oo) = 0.96(3). A third 
point is obtained from the effective mass a(l) = 1.09(2) 
in the fully paired superfluid phase ()19ej) . We use these to 
provide a smooth interpolation via the parametrization 

, hi . . 

z c = tanh — (24) 

as shown in Fig. [3] 

We use the variable z c (|2~4"]) here so that the interpo- 
lation is over a finite range z c 6 [—1,1]. The parameter 
c gives us some control over the shape of the resulting 
curve as demonstrated in Fig. |4j Note that for c < 2, 
dz c /dx\ x=Q = 0, hence the function a(x) will be flat at 
x = Q. The addition of a few spin-down particles should 
not affect the mass of the spin-up particle: most of the 
change should occur when one approaches equal densi- 
ties. There is not really sufficient information to further 
characterize this parametrization, but the effective mass 
does not vary much, so we do not expect this to be a 
significant source of error. 

Including density dependent inverse masses 
a a, bin a , nf,) can be of quantitative importance, but 
does not significantly alter the qualitative aspects of 
the ASLDA (such as presented in [lCj). Thus, to obtain 
qualitative results, it can be a good first approximation 
to simply use the functional with a a ^(n a , rib) = 1- 

D. "Hartree" Energy : f3(x). 



lation for the function G(x) (|23|) rather than directly for 
(3{x) or g(x) because: 1) It is finite everywhere and, 2) 
such an interpolation is independent of the inverse effec- 
tive mass function a(x). Again, we use the parametriza- 
tion ([24]) with z c 6 [—1,1] with the assumption that G(x) 
is smooth at x = 1, which requires that the interpolated 
function G(z) — G{x{z)) is smooth and even. This is 
quite easy to do (one can always just explicitly make 
the interpolation even \G(z) + G{— z)]/2). The assump- 
tion of smoothness also gives a non-trivial constraint on 
the data which should admit a smooth interpolation at 
z c = 0. The resulting function is shown in Fig. [5] for sev- 
eral values of c. Here the value of c = 2 is best because 
dz c /dx\ x= o = 2 is finite, preserving the structure of the 
interpolation. (The other values of c drastically affect 
the slope of the interpolated g(x) without any physical 
motivation.) The preference for c = 2 is most evident in 
Fig. El 




0.3 L"- 
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To finish the parametrization, we must provide an in- 
terpolation of the function g(x). We provide an interpo- 



1.1 



c 
II 

£ 1.0 



1 1 1 


c= 

c= 2 


^^^^ ^ 

. *N 

"'• 


N 

S 

^^^^ V 


1 1 1 1 1 l 



0.0 0.2 0.4 0.6 0.8 1.0 

x (lower curves), 1/x (upper curves) 



FIG. 4: (Color online) Effective mass as a function of x — 
rib/ria for several values of c: c = 1/2 - dashed green, c = 1 
- solid blue, c = 2 - dotted red. The upper curves represent 
the extension of the curves to x > 1 with abscissa 1/x. Since 
the effective mass is not likely to change rapidly near x — 
(especially for the majority component), the parameter c is 
probably best kept near 1 and we use c = 1. 



FIG. 5: (Color online) The function G(z) defined in ([23]) plot- 
ted for various values of c used to define the interpolation 
through the coordinate z c (|24[) : c = 1 - dashed green, c = 2 
- solid blue. See Fig. |S] which better emphasizes how this 
value is a better fit. We omitted the red points from our 
interpolation as discussed in the text. 



IV. LOFF 

Finally we include the derivation of a simple relation- 
ship between the average pressure and energy density of 
a LOFF state with period L imposed by boundary con- 
ditions. The physical solution will have a spontaneously 
determined length scale L, but our calculation must im- 
pose this through boundary conditions. To model these 
states, we use a periodic basis, summing over transverse 
and Bloch momenta. One must then adjust L to maxi- 
mize the average pressure. 

Consider the form of the average pressure at unitarity. 
The lack of scales implies 

s 5 V(s- 2 fi a , s- 2 fx b , sL) = T{na, fi b , L). (25) 
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TABLE I: Interpolation points for the function a(z) 



0.4 O.i 

X — Ub/Ua 



FIG. 6: (Color online) The function f(x) = g 5/3 (x). Both 
function plotted for various values of c used to define the 
interpolation through the coordinate z c ()24[) : c = 1 - dashed 
green, c = 2 - solid blue. It can be clearly seen here that the 
value of c = 2 is preferred for a smooth interpolation. 



1.000(00) 
1.094(17) 
0.962(28) 



TABLE II: Interpolation points for the function G(z) 



G 



-1.0000 
-0.6502 
-0.5886 
0.0000 
0.5886 
0.6502 
1.0000 



1.0000(00) 
0.6293(15) 
0.5797(20) 
0.3577(15) 
0.5797(20) 
0.6293(15) 
1.0000(00) 



Differentiating and using the thermodynamic relation- 
ship n a ,b = dV/dfi a ,b, we obtain: 



dV 

5V - 2^ a n a - 2 f i b n b + L— = 0. (26) 
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Coupled with the thermodynamic relationship V 
Ha.n a + Hbnb — £ we have 



dV 

X = L-V = 28- 3V. 
oL 



(27) 



When the pressure is maximized, this quantity X = 0. 
Note that this relationship is derived solely on dimen- 
sional grounds. 



V. INTERPOLATIONS: 



The functional (P4"a|) is completely described by the 
three functions at(x), f3(x), and j(x) through equa- 
tions ([8]). These functions are defined by the functions 
a(z) and G(z) that interpolate the data in tables U and 
ILT1 The resulting weighted cubic-spline interpolations are 
shown in Fig. [3] and Fig. [5]) respectively. 

a(x) = a (tanhlnx) , 

P{x) = 2 2 ' 3 G (tanh^) + 

a(x) a(x _1 ) 
~ (1 + xf/ 3 ~ (I + X- 1 ) 5 / 3 ' 

■y(x) = 7. 
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